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1 . SUMMARY 

A method has been developed which calculates two-dimen- 
sional, transonic, viscous flow in ducts. The finite volume, 
time marching formulation is used to obtain steady flow solu- 
tions of the Reynolds-averaged form of the Navier Stokes 
equations. The entire calculation is performed in the physi- 
cal domain. This paper investigates the introduction of a 
new formulation of the energy equation which gives improved 
transient behavior as the calculation converges. The effect 
of variable Prandtl number on the total temperature distribu- 
tion through the boundary layer is also investigated. 

A turbulent boundary layer in an adverse pressure gradi- 
ent (M ■ 0.55) is used to demonstrate the improved transient 
temperature distribution obtained when the new formulation of 
the energy equation is used. A flat plate turbulent boundary 
layer with a supersonic freestream Mach number of 2.8 is used 
to investigate the effect of Prandtl number on the dis- 
tribution of properties through the boundary layer. The 
computed total temperature distribution and recovery factor 
agree well with the measurements when a variable Prandtl 
number is used through the boundary layer. 
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2. INTRODUCTION 


This paper is an extension of the work reported else- 
where in this conference [1], A review of the features of 
the new method will be included here but a more complete 
discussion may be found in references 1 and 2. 

The features of the current method can be summarized as 
follows. Control volumes are chosen so that smoothing of 
flow properties, typically required for stability, is not 
needed. Different time steps are used in the different gov- 
erning equations to improve the convergence speed of the 
viscous calculations. A multi-volume method for pressure 
changes in the boundary layer allows calculations which use 
very long and thin control volumes (length/height ■ 1000). 

3. GOVERNING EQUATIONS 

The unsteady forms of the continuity equation, the x- 
momentum equation, the y-momentum equation, and the energy 
equation, in integral form, are used to obtain steady-state 
solutions for flow through 2-dimensional ducts. This ap- 
proach differs from our previous work [1] where the assump- 
tion of constant total temperature was used instead of the 
full energy equation. The ideal gas equation of state and a 
Prandtl mixing length turbulence model [1] complete the 
governing equations needed to solve for the unknown vari- 
ables p,u,v,P,p, and T. 

For a finite control volume where we can assign one 
value of density to the control volume, and for a finite time 
step, <St, continuity states that, 

P° + * - p° » Sp = -[ jj p u • d a] (1) 

where the integral is evaluated explicitly at the current 
time step, n. In arriving at an expression which relates the 
pressure change directly to the continuity error, we will 
assume that changes in temperature are small in comparison to 
other changes for one time step. Thus, we can relate changes 
in pressure to changes in density through the ideal gas equa- 
tion of state. 

p n+l _ P n = 6P = -RT[ jj p u • d A] (2) 

For the method introduced in the current work, a non-conserv- 
ative form of the unsteady momentum equation is used. The 
non-conservative form is used because it allows the current 
method to use different time steps for the continuity, momen- 
tum, and energy equations. The difference between the non- 
conservative and conservative forms of the unsteady momentum 
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equation is associated with the unsteady and convective 
terras. Specifically, we note that 

3 (pii) 3u 

" ■ g" t ■ “ + 7 • p 2.““ p Tt + p — * 7 u (3) 

and the right hand side of Eq. 3 can be rewritten as 
3u 3u 

p ^ + P u • V u “P jI + V • p u u - u (V • pu) (4) 

When the right hand side of Eq. 4 is combined with the pres- 
sure and viscous terms, the momentum equation in integral 
form becomes 

(u) n+1 - <u) n * 6(u) = [-// p u u • d A + u // p u • d A 


(5) 

- // Pfi^ • d A + // <y V u + u V u T ) ♦ d A ] 

To maintain stability, the properties must be updated in the 
proper sequence. In the current method, the sequence is: 

1. update the pressure from the continuity equation; 

2. update the velocities from the momentum equations using 
the new pressure and old velocities and old density; 

3. update the density from the ideal gas equation of state; 

4. update the temperature from the energy equation. 

4. ENERGY EQUATION 

For many calculations of transonic viscous flow, the 
assumption of constant total temperature will give a suffi- 
cient representation of the energy equation in the flow 
field. By assuming constant total temperature, the computa- 
tions are less expensive to run and the computer storage 
requirements are less. The assumption of constant total 
temperature is usually satisfactory if: 

1. an adiabatic wall is assumed in the calculations; 

2. no work Is done on the fluid at the solid boundaries; 

3. the Mach numbers in the flow fields are low enough that 
total temperature gradients within the boundary layer are 
small; 

4. the Prandtl number is approximately 1.0. 

For a Prandtl number of 0.9, the solution should not 
deviate greatly from the constant total temperature assump- 
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tion. However for high speed flow, the energy equation 
should be Included In the calculations especially If the 
Prandtl number deviates greatly from 1. 


Two forms of the Integral formulation of the energy 
equation will be derived next. 


The energy equation in differential form is 

3E T 

-y- + 7 • E^ u ™ -7 • + 7 • [u • (y 7 u + w V ^ 


where the total energy per unit volume, E t , is 


P (e + V 2 (u 2 + v 2 )) = pe t 


The left hand, side of Eq. 6 can be rewritten as 
3E 3 (pe ) 

TT + V * E t ii Tt— + 7 * p 


and 


3(pe t ) 

Tt 


+ V • (pe fc ) £ 


9e t 

TT 


+ p u • Ve 


P u 

( 6 ) 

(7) 

( 8 ) 

(9) 


then, expanding the right hand side of Eq. 9, we get, 

3e t 3e t 

p yy- + p u • 7e fc = p yy- + 7 u e fc - e fc (7 • p ji) (10) 

The procedure just outlined is identical to what was 
done to the unsteady and convective terms in the momentum 
equation (see Eqs. 3,4). 


The heat flux vector, can be represented as 

q = -kVT 


Substituting Eqs. 8-11 into Eq. 6 , we get 

3e 

p y—— * -7 • p u e fc + e^(7 • p u) - 7 


+ 7 


[ u • ( 11 V u + y 7 u ) ] - 7 


(-k7T) 
P u 


(ID 


( 12 ) 


The integral form of the energy equation is then 

5e. 


't 

TT 


p x 6 V 0 I = 


- // p u e t 


A + e / / P u • d A - // -kVT • d A 


(13) 
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+ // [ u • (|17 u + pV u T ) • d A - // P u • d A 

where e^, is an average value for the control volume. As 
with _ the momentum equation, Eq. 13 has a 
term e^ // p jj • d A , which removes the continuity error 
contribution To the” energy error. 

This form of the energy equation, when incorporated into 
the current method, behaved poorly. Initially there were 
large errors in continuity and momentum and these large er- 
rors acted through this energy equation to cause errors in 
the total energy for a control volume. This Interaction was 
destabilizing. 

An alternative form of the energy equation will now be 
derived. This alternative form has enhanced convergence 
properties when compared with the above formulation. Brief- 
ly, the energy equation is reformulated so that changes in 
total enthalpy, h t , are calculated rather than changes in 
total energy, e^, which was done previously. This allows us 
to see the terms which cause departures from uniform total 
temperature - for both the steady state solution and the 
transient solution. 


The total enthalpy can be defined in terms of the total 
energy and the static temperature 


or 


e fc + P/p 


+ RT 


(14) 

(15) 


Taking the derivative with respect to time and multiplying by 
the density, we get 

3h^ 3e. 

r wr 

(16) 


t 

p Tt- 


3T 


p TT + p R Ft 


The static temperature T can be represented in terms of the 
total enthalpy and the absolute velocity as 


Therefore 



t V 

(17) 


T a — 

<r 2c 

p p 

3T 

i 3h t _ V_ 3V 

(18) 

3 1 

c at c at 

p p 

18 into 

Eq. 16, we obtain 


3e t 

0 Sh t R „ 3V 

(19) 

p TT 

= Y TT + 0 TT v Tt 

P 
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where y is the ratio of specific heat capacities and V is the 
magnitude of the velocity vector. Using equations (19) and 
(14) to eliminate e fc from equation (12) we get 
3h 

| * -7«p u h fc + (h fc - £-)(V*p u) + 7‘kVT 

( 20 ) 

+ 7*[u • (p7u + u7 u T ] - p V || 

P 

Using h t ■ C T + V^/2 and k “ u C /Pr, k7T may be replaced 
by F p 

2 

k7 T » Vh t - ^ 7 (£-) (21) 

and from continuity we may replace 7*p u with -3p/3t. 

Therefore the energy equation written as a conservation equa- 
tion for total enthalpy is 
3h 

- = " 7 * p “ h t + h t ^ 7 * p “) + 7 *^ 7h t 

(I) (II) (III) (22) 

+ 7 -"(i - + 7 -*‘<r 7 > ; ♦ IS? - £ 7 H 

P 

(IV) (V) (VI) (VII) 

Terms I and II when combined give - p u • Vh t . Therefore 
terms I + II and III contain h fc only in tRe form Vh . Thus, 
when these are the only important terms in the equation, flow 
with uniform total temperature at the inlet will retain this 
uniform total temperature provided that the boundary condi- 
tions are consistent with this. 

Term IV is a viscous transport term for total enthalpy when 
the Prandtl number is other than 1. Term V is another vis- 
cous transport term. It however contains the expression 
(u • 7) u which is the gradient of the velocity in the direc- 
tion of the velocity; these gradients are usually small com- 
pared with other velocity gradients. Since terms IV and V 
have the form 7 • ( ), they are not source terms, rather they 
can only redistribute the total enthalpy. Terms VI and VII 
on the other hand have the form of source terms. Relative to 
the steady state, they are proportional to the continuity 
error and the momentum error, respectively. We may write 
them as 


M 


l 


3p 

Tt 


, 3V 

+ m Tt 


At the steady state, Eq. 22 becomes 


(23) 
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0 


( 24 ) 


2 

« -V’p u h t + 7 • 7h fc + V*y(l - ^)V(— ) 

+ V • y(u • V) u 

Therefore we may artibrarily alter the variables 1 and m in 
Eq. 23 and the steady form of the energy equation, Eq. 24, 
will be obtained for converged solutions. The transient 
behavior of h t is improved in the calculation procedure by 
choosing 1 * m * 0, i.e. by omitting the transient source 
terms in the enthalpy equation. 

In integral form then the equation for enthalpy changes is 
6h t 

P -£■£— 6Vol =■ y{- // p u h fc • d A + h // p u • d A 


(25) 

+ // Vh t •. d A + // [(u - iL) u»7u T + yu‘Vu]*dA 


where u a u„ + 


M t and w 


*5 


TT 


The time step used for the enthalpy equation is the same as 
for the momentum equation. If the transient source 
P 3p 

term — had been retained in the enthalpy equation, it 

would have been necessary to link the continuity and energy 
equation time steps. Omitting this term allows us to use 
different time steps for the energy equation. 


5. TEST CASES 


Two test cases will be used to explore various aspects 
of the more complete form of the energy equation, Eq. 25, 
discussed previously. 

5.1 Turbulent Boundary Layer in an Adverse Pressure Gradient 

The geometry and grid used in this test case are shown 
in Fig. 1. Flow in this geometry was used in Ref. 1 to test 
the accuracy of the new computational scheme. In Ref. 1 the 
velocities in the duct were low enough that the flow could be 
treated as incompressible. Here, the inlet freestream Mach 
number was increased to 0.55. The purpose of this test case 
was to illustrate the advantage of the new formulation of the 
energy equation. 

The static temperatures presented in Fig. 2 are from 
calculations after 500 iterations. It can be clearly seen 
that the new formulation, Eq. 25, gives a better transient 
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solution to the energy equation and it should result in a 
reduction in the computer time required to reach a steady 
state solution. Fig. 3 shows the corresponding total temper- 
ature profiles for the two formulations of the energy equa- 
tion. 



Fig. 1 Grid and Geometry Used to Demonstrate the 
Advantages of the New Formulation of the 
Energy Equation 



Fig. 2 Static Temperature Distribution Through 
the Boundary Layer at M=0.55, x=200 mm, 
after 500 iterations 

5.2 Flat Plate Turbulent Boundary Layer at M * 2.8 

Van Driest [3] presents the total temperature distribu- 
tion within a flat plate turbulent boundary layer with a 
freestream Mach number of 2.8. The experimental total tem- 
perature distribution is shown in Fig. 4. The geometry and 
grid for these calculations are shown in Fig. 5. The height 
of the duct was 63.5 mm and the length of the duct was 254 
mm. The computational grid shown in Fig. 5 consists of 21 
axial grid points and 14 transverse grid points. The inlet 
boundary layer thickness of 6.35 ram was 10% of the duct 
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Fig. 3 Total Temperature Distribution Through the 
Boundary Layer at M=0.55, x=200 mm, 
after 500 iterations 

height. The Reynolds number based upon axial distance is 
approximately 10'. To stabilize these supersonic flow 
calculations, the upwind effective density method was used 
[2]. This means that an effective density used at a grid 
point is calculated with the ideal gas equation of state 
using the pressure from the next upstream grid point. The 
inlet velocity, total temperature, and total pressure were 
specified at the upstream boundary. Three calculations were 
performed with different assumptions about the turbulent 
Prandtl number. These assumptions were 

1. Pr t = 0.90 Pr. - 0.73 

2. Pr t - 0.73 Pr* = 0.73 

3. Pr t varies linearly through the boundary layer from 0.9 
at the wall to 0.66 in the freestream. 


The turbulent Prandtl number is typically set equal to a 
constant of 0.9 in calculations [4]. The calculated total 
temperature distribution through this boundary layer using a 
constant turbulent Prandtl number of 0.9 is shown in Fig. 6 
(represented as Q ) . The recovery factor is calcu- 

lated to be 0.920 which compares with the empirically deter- 
mined value of 0.90. However, the distribution of total 
temperature through the boundary layer does not compare well 
with the experiment. If the turbulent Prandtl number is set 
equal to the laminar Prandtl number of 0.73, the total tem- 
perature distribution changes as seen in Fig. 6 (represented 
at O's). The distribution through the outer part of the 
boundary layer has improved but the recovery factor of 0.813 
does not compare well with the experimental value of 0.90. 
Schlichting [5] notes that the turbulent Prandtl number is 
not constant through the boundary layer. The experiments of 
H. Ludwieg [6] for turbulent flow through a pipe show that 
the Prandtl number varies from approximately 0.9 at the pipe 
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Fig. 4 Experimental Total Temperature Distribution 
in a Flat Plate Turbulent Boundary Layer 
M=2.8, Re = 10? after van Driest 



Fig. 5 Geometry and Grid For Boundary Layer Calculations 
at M=2.8 



Fig. 6 Total Temperature Distribution For Flat Plate 
Boundary Layer at M=2.8 




wall to 0.66 at the center of the pipe. This distribution Is 
shown In Fig. 7. The variation is almost linear. For the 
third set of calculations, the Prandtl number was assumed to 
vary linearly through the boundary layer from 0.9 at the wall 
to 0.66 at the edge of the boundary layer. The total temper- 
ature distribution for this case Is shown In Fig. 6 (repre- 
sented asA's). The total temperature distribution calcu- 
lated using a variable Prandtl number is also compared with 
the experimental results in Fig. 8. Both the distribution of 
total temperature through the boundary layer and the recovery 
factor of 0.90 are in good agreement with the experimentally 
measured values. 



Fig. 7 Ratio of the Turbulent Transfer Coefficient Over 
the Length of a Radius in Turbulent Pipe Flow 



Fig. 8 Total Temperature Distribution For Flat Plate 

Boundary Layer at M*2.8 Computation vs. Experiment 
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6. CONCLUSIONS 


A new formulation for the energy equation was introduced 
which has improved transient behavior when compared with the 
standard formulation. The new formulation removes the in- 
fluences of continuity and momentum errors from the energy 
equation during transients in the solution. 

For flat plate turbulent boundary layer flow with a 
freestream Mach number of 2.8, the calculated total tempera- 
ture profile was improved by using a variable Prandtl number 
through boundary layer. The recovery factor of 0.90 agreed 
very well with the empirically determined value of 0.9. 
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